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1.  The  density  matrix  from  a  numerical  solution  of  (3)  with  the  Hamiltonian 

given  by  (22)  and  (31),  the  dissipation  is  given  by  (23) . 

2.  The  density  matrix  as  in  Figure  1  but  with  y  =  1/100 . 


1.  Introduction 


It  is  common  that  the  differential  equations  used  to  model  physical  systems  possess 
some  “structural”  features  that  embody  important  physical  properties  (e.g.,  differential- 
algebraic  invariants  such  as  energy)  which  are  reflected  in  the  qualitative  behavior 
of  the  system.  Often,  this  structure  is  purposely  built-in  to  the  model  on  physical 
grounds.  Unless  the  numerical  methods  are  specially  designed,  a  numerical  solution  of 
such  a  system  will  not,  in  general,  respect  this  structure,  leaving  open  the  possibility 
of  qualitatively  incorrect  results.  Structural  considerations  are  particularly  important 
when  the  time  domain  of  interest  is  much  larger  than  the  system’s  characteristic  time 
scale(s).  In  such  cases  it  can  be  difficult,  if  not  impossible,  to  obtain  even  qualitatively 
correct  results  when  the  numerical  method  does  not  preserve  the  system’s  structure.  (A 
particularly  interesting  example  is  the  simulation  of  the  long-time  evolution  of  the  solar 
system  where  phase-space  conserving  (symplectic)  techniques  proved  indispensable  [1].) 
At  the  very  least,  the  penalty  for  ignoring  the  essential  structure  of  the  system  is  a 
greatly  increased  computational  cost  forced  by  the  necessity  of  a  small  time  step  to 
keep  the  effects  of  structural  errors  from  accumulating  catastrophically.  For  a  more 
detailed  discussion  of  these  effects  along  with  various  examples  see  Shadwick  et  al.  [2] 
and  the  references  therein. 

The  goal  is  not  to  eliminate  all  numerical  error  but  rather  to  identify  and  minimize 
those  classes  of  numerical  errors  that  are  most  damaging  to  the  solution.  Backwards 
error  analysis  [3-7]  provides  significant  understanding  in  how  a  numerical  method  may 
lead  to  a  numerical  solution  with  qualitatively  incorrect  behaviour.  The  basic  approach 
is  to  take  the  view  that  when  a  numerical  method  is  applied  to  a  given  differential 
equation,  the  resulting  numerical  solution,  while  an  approximation  to  the  exact  solution 
of  the  original  equation,  is,  in  fact,  the  exact  solution  of  some  differential  equation. 
Clearly  this  second  differential  equation  is  related  to  the  original  equation,  and  the 
form  of  this  new  equation  can  often  lead  to  valuable  insight  into  the  behaviour  of  the 
numerical  method  [8].  Of  particular  importance  will  be  the  physical  consequence  of  any 
coupling  terms  not  present  in  the  original  system;  it  is  terms  of  this  type  which  are  the 
cause  of  structural  errors. 

These  ideas  are  also  relevant  to  weakly  non-ideal  systems,  for  example  systems 
with  a  small  amount  of  dissipation.  Using  the  technique  of  operator  splitting  [9, 10], 
one  separates  the  differential  operator  into  two  pieces:  the  operator  representing  the 
ideal  part  of  the  system  and  whatever  non-ideal  terms  remain.  A  structure-preserving 
integrator  is  used  for  the  former  and  a  generic  method  for  the  latter.  This  technique 
guarantees  that  the  numerical  solution  has  the  proper  limit  as  the  dissipation  terms 
vanish.  Specifically,  consider  the  system  of  equations 

ij)  =  City]  +  C2[ (A]  (1) 

and  let  <Si(t)  and  ^(r)  be  the  respective  numerical  evolution  operators  for  a  time 
interval  r.  If  Si  and  <S2  are  accurate  to  at  least  second  order,  then  a  second  order 
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(2) 


approximation  to  the  full  evolution  is 

i>(t  +  t)  =  <Si(|  r)  o <S2(r)  o  r)ip(t)  +  0(t3)  , 

where  o  denotes  operator  composition. 

Consider  a  weakly  dissipative  system,  which  in  the  absence  of  dissipation  possesses 
some  structure  (for  example,  in  the  ideal  limit  the  system  may  be  Hamiltonian).  In 
this  case,  the  structural  errors  committed  by  standard  methods  can  mimic  the  effects 
of  dissipation,  enhancing  (or  perhaps  masking)  what  would  otherwise  be  a  small  effect. 
In  strongly  non-ideal  systems,  the  errors  introduced  in  evolving  the  ideal  part  of  the 
system  are  not  as  significant  in  that  they  are  not  likely  to  yield  behaviours  that  are 
not  otherwise  present.  The  same  is  true  for  errors  encountered  in  the  evolution  of  the 
dissipative  terms.  Empirically,  one  finds  that  dissipation  terms  tend  to  be  sufficiently 
simple  that  they  possess  no  delicate  structure.  Consequently,  numerical  errors  will 
generally  not  result  in  qualitatively  different  behaviour,  merely  in  different  effective 
values  for  the  dissipation  parameters. 

2.  n-Level  Quantum  Systems  with  Weak  Dissipation 

Our  interest  here  is  to  explore  how  these  ideas  can  be  applied  to  an  n-level  quantum 
system  described  by  a  density  matrix,  p,  subject  to  both  reversible  and  irreversible 
processes.  A  general  dynamical  equation  governing  the  evolution  of  p  is  the  master 
equation  [11] 

hp=-i  {H,p]+A{p\,  (3) 

where  the  hermitian  Hamiltonian  H  describes  the  reversible  dynamics  and  the  (general) 
dissipation+  is  represented  by  the  linear  operator  A.  Examples  include  an  atomic 
system  interacting  with  an  external  radiation  field  subject  to  spontaneous  emission  or  a 
collection  of  ions  in  a  trapping  potential  subject  to  external  laser  pulses  and  decoherence. 
We  are  most  concerned  with  the  case  where  the  irreversible  processes  are  weak  compared 
to  the  reversible  processes,  i.e.,  the  case  ||A||  -C  ||ffj|. 

For  the  moment,  let  us  consider  the  properties  of  the  quantum  Liouville 
equation  [12],  the  ideal  counterpart  of  (3): 

hp=-i[H,p\.  (4) 

This  equation  has  a  non-trivial  kinematic  structure  —  the  Hioe-Eberly  invariants 

Cj  =  tr  (P ,  j  =  1, 2 . . .  n  (5) 

are  non-evolving  regardless  of  the  form  of  the  Hamiltonian  [13].  These  invariants  are  a 
direct  consequence  of  the  unitary  evolution  of  the  density  matrix  and  are  the  analogues 
of  the  Poincare  invariants  in  classical  mechanics,  carrying  information  of  equal  physical 

+  Although  here  we  are  primarily  concerned  with  dissipative  processes,  A  can  represent  arbitrary  non- 
Hamiltonian  processes. 
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import.  A  numerical  solution  of  (4)  where  the  Cj  are  not  preserved  is  thus  in  danger  of 
being  unphysical. 

In  the  method  of  unitary  integration  [14],  the  structure  of  the  quantum  Liouville 
equation  is  preserved  exactly  since  the  numerical  time  advance  map  is  constructed  as  a 
unitary  transformation.  That  is,  the  density  matrix  is  evolved  through  a  time  interval 
r  according  to 

Qn+l  =  UQnU\  (6) 

where  gn  is  the  numerical  approximation  to  p(n  r )  and  U  is  a  unitary  matrix  which 
approximates  the  true  evolution  operator  4-  t)  accurate  through  order  k  in  r: 

U  =  il(t,  t  +  r)  +  0(tk+1)  .  (7) 

Thus,  although  gn  is  only  a  numerical  approximation  to  the  true  solution  of  the  quantum 
Liouville  equation,  the  Cj  are  exactly  conserved.  Furthermore,  since  the  numerical 
method  advances  gn  in  time  by  a  unitary  transformation,  we  are  guaranteed  that  the 
numerical  solution  will  be  hermitian. 

Returning  to  (3),  our  strategy  is  to  use  a  unitary  integrator  for  the  Hamiltonian 
part  of  the  evolution  and  a  generic  algorithm  for  the  dissipative  terms.  Thus  in  the 
limit  ||A||  — ►  0,  the  Hioe-Eberly  invariants  are  exactly  conserved.  In  this  way  we 
guarantee  that  all  dissipation  and  decoherence  (variation  of  the  Cj)  are  due  to  the  non- 
Hamiltonian  terms,  and  not  to  any  numerical  artifacts. 

3.  Constructing  Unitary  Integrators 

As  with  every  unitary  matrix,  U  can  be  expressed  as  the  exponential  of  an  anti-hermitian 
matrix.  For  our  purposes,  it  is  convenient  to  write 

U(t,r)  =  e-iTA{t'T),  (8) 

where  A  is  hermitian.  The  matrix  A  is  computed  by  expanding  (6)  in  a  Taylor  series 
about  r  —  0  and  matching  term-by-term  with  the  Taylor  series  for  p(t  +  r)  obtained 
from  the  quantum  Liouville  equation.  A  straightforward,  if  tedious,  calculation  [14] 
reveals  the  following  approximations  for  A:  to  second  order 

A  =  H(t)  +  ^rH'(t)-  (9) 

to  third  order 

A  =  H(t)  +  i  r  H'(t)  +  i  r2  H"(t)  +  ^  r2  [H(t),  ;  (10) 

and  to  fourth  order* 

4 = h  (») + if  h  ((>' + ir  2h  M" + if sm" 

+  ~r‘  [ H(t ),  Hftf]  +  if3  [tf  (f),  H(tf} .  (11) 

*  There  is  a  misprint  in  (6)  of  Ref.  [14]  which  is  corrected  here. 
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Note  that  to  obtain  accuracy  beyond  second  order,  one  must  take  into  account  that,  in 
general,  [Hfa),Hfa)]  ±  0. 

Fortunately,  to  use  these  expressions  for  U,  it  is  not  necessary  to  exponentiate  a 
(general,  possibly  quite  dense)  n  x  n  matrix.  We  are  free  to  approximate  A  in  any 
way  consistent  with  the  order  of  the  method.  This  freedom  can  be  used  to  simplify  the 
construction  of  U(t,  r).  In  the  following,  it  is  convenient  to  assume  that  H  is  traceless. 
To  begin,  let  \k  be  a  basis  for  the  set  of  all  (traceless)  n  x  n  hermitian  matrices,  with 
normalization 

tr  A*  Xj  =  2  8ij,  (12) 


whence 


n2- 1 

A  =  ^  ak  A*, ,  where  ak  =  \  tr  A  Xk .  (13) 

k=l 

In  addition,  we  may  choose  the  A*  in  such  a  way  that  either  Xk  is  diagonal  or  has  at 
most  two  non-zero  elements.  In  either  case,  e'aXk  is  easily  computed.  Our  goal  is  to 
replace  the  single  exponential  in  (8)  with  a  product  of  exponentials  of  the  A*: 


n2- 1 

e~iTA  =  e-iT^Afc  +  0(tk+1),  (14) 

k=i 

for  some  (3k.  Since  the  Xk  do  not  necessarily  commute  with  one  another,  the  j3k  will 
depend  on  the  aj  in  a  complicated  manner. 

While  it  is  possible  to  determine  the  /3k  in  (14)  by  repeated  application  of  either 
the  Campbell-Baker-Hausdorf  [15-17]  or  the  Zassenhaus  [18]  formula,  such  an  approach 
quickly  becomes  quite  complicated  as  n  increases  and  does  not  seem  to  be  well-suited 
to  symbolic  computation.  A  more  straightforward  approach  is,  using  (13),  to  expand 
the  left  and  right  hand  sides  of  (14)  as  polynomials  in  r  and  match  terms  order  by 
order  [19].  For  an  integrator  of  order  /c,  we  can  write  ak  as 

ak  =  '%2Ti<xk)  >  (15) 

l=o 

and  similarly 

A  =  £^  +  0(0-  (16) 

j=0 

Using  these  expressions  in  (14),  we  obtain 

1  -  i  t  rj  dp  Xk  -  \  r2  53  2Z  Tj+j’ak)  ak']  xk  Afc' 

1=0  k= 1  1,1'=0  k,k'=l 

0® ft* ’ Aj  + 

k=\  \  1=0  1,1'= 0 
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Collecting  powers  of  r  we  find 

-  It  £  of  A*  -  ir’  £  of*.  -  Jt*  £  (a<">)2  A; 


n2-l 


fc=l 


*=1 


k=l 
n2- 1 


~  5 1-2  5Z  ai0)  al°}  Afc  Afc/  +  . . . 

k^k' 

II  f  1  -  i^°’  ^  -  ir2  >  A,  -  1  r2  (/S®)’  A?  +  ...)-  1 

fc=l  '  ' 

=  -ir  §/£  (tf')2  A| 

ft=l 


fc=l 


k~\ 


(18) 


k‘>k 


Multiplying  (18)  by  Ap  and  using  (12)  yields 


n2- 1 


r  a£0)  +  T2a^  -  j  r2  ^  tr  AP(A*  Afc<  +  A*<  Afe)  +  . . . 

k’>k 

=  r  $°>  +  r2/?W  -  i  r2  ]T  /?f  /#>  tr  Ap  A,  V  +  •  •  •  • 


k'>k 


Equating  like  powers  of  r  we  have 

/5f  =  af; 


.  n2— 1 


J  ak]  °i?  tr  AP  tA*>  Afc'l ; 


(19) 

(20) 
(21) 


k’>k 


/K_  i) 

Clearly,  one  can  continue  this  expansion  to  obtain  equations  for  frp  ...  /?p  .  Moreover 

this  technique  is  well  suited  to  straightforward  implementation  in  a  symbolic  algebra 
language.  In  such  an  implementation,  one  uses  the  explicit  form  of  the  A*,  to  obtain  a 
matrix  equation  at  each  order,  which  naturally  leads  to  a  set  of  linear  equations  for  p[j\ 
This  approach  eliminates  the  need  for  developing  explicit  expressions  for  {}jf\ 


4.  An  Example 


By  way  of  example,  we  consider  a  two-level  atom  with  population-preserving  decay  [11]. 
The  Hamiltonian  for  this  system  is 


H  = 


Q(t) 


n(t)  ' 

— e 


(22) 
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where,  for  simplicity,  we  have  taken  0(f)  to  be  real.  The  dissipation  operator  is  given 

by 


A  \p]  = 


1P22 


-5TP12  1 

-1P22  J  ’ 


(23) 


As  we  explained  above,  we  are  interested  in  the  case  of  weak  dissipation,  i.e.,  in  7  <C  e. 

Following  the  notation  of  (13),  for  a  second  order  integrator  we  have  ch=Q.+\t  O', 
oi2  =  0  and  a3  =  e,  and  we  need  to  compute  ft  through  order  r.  From  (20)  and  (21) 
we  find  ft  =  Q  4-  \  r  O',  ft  =  r  e  fi,  and  ft  =  e.  Explicitly  the  unitary  integrator  is 

r  cos(Qr  +  |fl'r2)  -i  sin(flr  +  \ Sl'r2)  " 

~~  k  -i  sin(Qr  +  |fl'r2)  cos(Or  +  |fl'r2)  ^ 


(  cos (efir2)  -sin(eQr2)  |  f  e  leT  0  1 
[  sin(eflr2)  cos(eQr2)  J  [  0  e,£T  J 

We  obtain  an  evolution  operator,  </(r),  for  the  dissipative  part  by  solving 


P  =  A[p] 


(25) 


using  the  midpoint  rule: 

ft,+i  -  en  =  t  A[|  (ft,+i  +  ft,)]  •  (26) 

We  can  interpret  (26)  as  differencing  (25)  at  t  +  \r,  which  yields  an  algorithm  accurate 
to  second  order.  Since  A  is  a  linear  operator,  (26)  can  be  solved  directly  to  provide  an 
explicit  formula  for  ft,+i): 


6n+l 


(fti)ll  + 


r  7 


1  +  |t7 

1  - \t  7 


(6n) 


22 


1  +  3T7 


(fc> 


!-3t7  ,  \  ' 

TT1T7  &  12 

1  -  §t7 


21 


1  +  5  r  7 


(ft.) 


22 


Following  (2),  we  can  write  the  total  evolution  operator  as 


where 


ft,+  l  =  T)  O  U(t,  T)  O  Q{\ t)  ft,  . 


W(t,r)  ft,  =  t/ft,^1 , 


(27) 


(28) 

(29) 


is  the  unitary  integrator.  Note  that  both  tr  Q  —  1  and  tr  U  =  1  and  hence  tr  ft,+i  =  tr  ft,. 

Below  we  compare  the  results  of  using  this  split-operator  method  with  the  results 
of  a  second-order  predictor-corrector  algorithm  (Heun’s  method  [20]).  Applying  the 
predictor-corrector  formulae  to  (3)  gives 

Qn+i  =  Qn-ir  [H (nr )  +  H((n  +  1  )r),  Qn]  +  r  A[ftJ  +  \  r2A  [A[ft,]] 

-  \ r2  [H((n  +  1  )r),  [H(nr),  ft,]]  -  \  ir2A^[H((n  +  1  )r),  ft,]] 

-|ir2[ff((n+l)r),A[e„]],  (30) 
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where  we  have  made  use  of  the  fact  that  A  is  linear.  Now,  as  expected,  tr  £n+i— tr  gn  =  0, 
since  the  trace  of  a  commutator  vanishes  as  does  tr  A,  thus  this  algorithm  also  exactly 
conserves  the  population. 


Figure  1.  The  density  matrix  from  a  numerical  solution  of  (3)  with  the  Hamiltonian 
given  by  (22)  and  (31),  the  dissipation  is  given  by  (23),  where  e  =  1,  Qo  =  *o  =  10, 
cr=^,7=^j,r=^j.  The  initial  condition  is  given  by  (33).  The  heavy  dashes  and 
solid  lines  denote  the  predictor-corrector  and  unitary  integrator  solutions  respectively. 
The  panels  are:  inversion  ((gn) 22  -  (Qn)  11)  (a),  error  in  the  inversion  (6),  Im  (gn)i2  (c), 
error  in  Im(^n)i2  (d),  tr^  (e),  error  in  trp£  (/),  energy  expectation  value  (tr  Hgn) 
( g ),  and  energy  error  ( h ). 
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Figure  2.  The  density  matrix  as  in  Figure  1  but  with  7=  ioo- 


We  provide  a  pair  of  numerical  examples  to  illustrate  the  advantage  of  the  structure- 
preserving  approach  over  the  general-purpose  algorithm.  We  take  the  interaction  to  be 
a  gaussian  pulse: 

a(i)  =  Q0e-(t-to)2/<72,  (31) 

and  use  the  parameters: 

e  =  1,  =  §,  <o  =  10,  o  =  f  and  7  =  ^- 


(32) 


We  set  the  initial  condition 


p(t  =  0)  = 


1 

2 

1  p— i7r/4 


I  pi7r/4 
4  e 


(33) 


and  use  the  same  time  step  r  =  ^  for  both  methods.  The  first  set  of  numerical  results 
is  shown  in  Figure  1.  Clearly  the  unitary  integrator  produces  smaller  errors  in  the 
inversion  and  coherence  elements,  and  is  also  somewhat  better  in  the  energy  than  pre¬ 
dictor-corrector.  As  with  symplectic  integrators,  the  nature  of  the  energy  error  is  not 
surprising.  Since  the  unitary  integrator  is  effectively  constructed  from  an  approximate 
Hamiltonian,  we  cannot  expect  tr  Hgn  to  be  computed  exactly  but  rather  we  expect  the 
error  in  tr  Hgn  to  oscillate  rather  than  exhibit  secular  growth  [14]. 

Figure  2  shows  numerical  results  with  the  same  parameters  as  above  except 
for  7  =  1/100.  In  both  cases,  the  energy  error  for  both  methods  is  quite  similar  and 
provides  little  insight  into  which  method  provides  the  more  accurate  solution.  The 


errors  in  tr  p2  are  more  telling:  the  operator  splitting  method  clearly  performs  better. 
In  Figure  2(e),  we  see  how  the  structural  errors  introduced  by  the  predictor-corrector 
algorithm  reduce  the  decoherence  caused  by  the  dissipation. 


4-1.  Discussion 


In  view  of  the  results  shown,  we  considered  several  different  approaches  to  handling 
the  evolution  of  the  dissipative  operator,  in  order  to  understand  the  origin  of  the 
superior  performance  of  the  algorithm  based  on  the  unitary  integrator.  In  addition 
to  using  the  midpoint  rule,  splittings  based  on  a  predictor-corrector  solution  to  (25)  as 
well  as  splittings  based  on  the  exact  solution  of  (25)  were  explored.  (Since  the  exact 
solution  requires  the  evaluation  of  an  exponential,  it  is  much  more  efficient  to  solve 
the  differential  equation  than  to  employ  the  exact  solution.)  In  each  of  these  cases  the 
full  evolution  was  computed  via  (28)  with  the  appropriate  choice  for  Q.  Regardless 
of  the  method  of  handling  the  dissipation,  essentially  the  same  results  were  obtained; 
clearly  the  advantage  of  the  operator  splitting  approach  is  due  to  the  use  of  the  unitary 
integrator. 

To  understand  these  results  better,  it  is  worthwhile  to  examine  the  numerical 
methods  in  detail  through  backwards  error  analysis  [3-7].  In  Section  1,  we  argued  that 
dissipation  operators  are  typically  sufficiently  simple  to  not  possess  a  “structure”  that 
is  physically  important.  This  assertion  is  borne  out  by  the  true  dynamical  equations 
corresponding  to  the  various  numerical  solutions.  For  the  full  system  (3),  this  analysis 
is  rather  cumbersome,  so  we  consider  separately  the  Hamiltonian  and  dissipative  terms. 

We  begin  by  examining  the  consequence  of  applying  the  midpoint  rule  to 


p  =  Aeff[p] 


[  7i  P22 

—  |  72^12 

[  -5  72P21 

-7i  P22  J 

(34) 


The  idea  is  to  determine  71  and  72  such  that  when  the  midpoint  rule  is  applied  to  (34) 
we  obtain  the  exact  solution  of  (25).  Since  the  midpoint  rule  is  second  order,  we 
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expect  71,2  =  7  +  0(r2).  The  procedure  is  to  compare  the  midpoint  rule  solution  of  (34) 
with  the  exact  solution  of  (25),  generated  from  the  Taylor  series  expansion  about  r  =  0. 
Some  straightforward  algebra  gives 

71  =  7  —  h  ^  +  lio  ^  +  ’  (35) 

12=1  ~h> 1-273 + 1^0  ^ + (36) 

What  this  says  is  that  when  we  use  the  midpoint  rule  to  solve  (25)  we  are  in  fact  getting 
the  exact  solution  of  the  slightly  modified  system  (34).  Now  although  71  ^  72,  the 
physical  meaning  of  the  dissipation  operator  is  only  slightly  changed  since 

21  =  1  — IrV  +  0(tV).  (37) 

72  16 

The  important  point  is  that  this  modified  dissipation  operator  is  physically  reasonable, 
with  qualitative  and  quantitative  behaviour  very  similar  to  the  original  dissipation 
mechanism.  This  is  precisely  what  it  means  not  to  have  a  delicate  structure;  the 
quantitative  and  qualitative  behaviour  is  robust  to  small  perturbations  in  the  form 
of  the  dynamical  equations.  However,  had  the  dissipation  model  been  more  complex, 
it  is  quite  likely  that  to  determine  Aeff  it  would  have  been  necessary  to  introduce  new 
couplings  in  the  equations  which  would  most  likely  lead  to  qualitatively  different  physical 
content.  This  is  exactly  what  we  will  see  below  when  we  perform  a  similar  analysis  on 
the  Hamiltonian  part  of  the  (3). 

Performing  this  same  analysis  for  the  predictor-corrector  solution  of  (25)  yields  a 
very  similar  result.  In  this  case  we  find 

71  =7  +  ^t273  +  ^t374  +  r4  7s  I  0(rr’y€),  (38) 

72  -  7  +  i  r V  +  i  T3  7<  +  i  r<  y  +  0(rs7») ,  (39) 


and 

—  =  1  -  l  rV  +  C>(r373) ,  (40) 

72  8 

leading  us  to  the  same  basic  conclusion  as  for  the  midpoint  algorithm. 

We  now  turn  to  an  analysis  of  the  numerical  solution  of  the  purely  Hamiltonian 
system  (4).  For  this  system,  the  predictor-corrector  time  advance  is 

Qn+ 1  =  Qn-iT  [H{nr)  +  H{(n  +  1)t),  ftj  -  \  t2  [H((n  +  1  )t),  [H(nr),  ft,]]  (41) 


In  the  same  fashion  as  above,  we  seek  to  determine  the  effective  Hamiltonian 


( 


^eff 


(42) 


such  that  (41)  applied  using  He$  gives  the  exact  solution  to  (4).  It  turns  out  that  for 
predictor— corrector,  it  is  possible  to  achieve  this  goal  only  through  third  order  in  t;  the 
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(43) 


fourth  order  corrections  result  in  a  non-hermitian  He^.  To  third  order  we  find: 


feff  =  e 
f2eff  = 


’i-|T»(n»  +  e>)] 
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o 


(44) 


At  fourth-order,  the  effective  system  is  no  longer  Hamiltonian  but  contains  dissipation, 
thus  we  expect  the  predictor-corrector  to  introduce  errors  in  the  Hioe-Eberly  invariants 
at  fourth-order  (for  example,  see  (11)  in  Ref.  [14]).  This  calculation  becomes  increasingly 
complex;  we  can  adequately  illustrate  our  point  more  clearly  by  considering  a  constant 
Hamiltonian.  At  fourth-order,  solving  (4)  with  (41)  is  equivalent  to  the  exact  solution 
of 


hp=  —  i  [Hefi,p]  +  Aeff[p] , 


(45) 


with 


and 


=  € 

fleff  ~  Cl 


l-p(fi’  +  62) 


+  i  r3e  Cl  (Cl2  -I-  e2)  , 


Aeff[p]  =  r3  Cl2 


P22  ~  Pll 
Pl2  ~  P21 


P2\  ~  Pl2 
Pll  —  P22 


(46) 

(47) 


+  2r3e 


f2  (Q2  +  e2 


r 

) 


0 

Pll  ~  P22 


Pll  ~  P22 
0 


—  2  r3  e2  (Q2  +  e2)  |  ^  j  -  (48) 

Notice  that  the  form  of  Aeff  describes  a  much  different  set  of  physical  processes 
than  those  responsible  for  the  actual  dissipation  in  (3).  For  example,  the  presence  of 
the  inversion  in  the  off-diagonal  elements  and  the  fact  that  Aeff  drives  the  system  towards 
equal  populations  rather  than  decay  to  the  lower  state.  Although  ||Aeff||  <£.  ||A||,  the 
numerical  results  show  that  Aeff  is  a  sufficient  perturbation  to  significantly  affect  the 
solution.  Thus  even  though  (3)  does  not  exactly  preserve  the  Cj,  the  remnants  of  the 
ideal  structure  of  (4)  are  sufficiently  strong  that  the  perturbations  introduced  by  the 
predictor-corrector  are  enough  to  adversely  affect  the  accuracy  of  the  numerical  solution. 


5.  Conclusions 

The  advantages  of  unitary  integration  are  seen  most  dramatically  in  ideal  or  nearly  ideal 
systems  sensitive  to  small  errors.  The  important  point  to  be  made  here,  however,  is  that 
the  operator  splitting  approach  where  a  unitary  integrator  is  used  for  the  Hamiltonian 
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evolution  is  guaranteed  to  recover  the  ideal  solution  in  the  limit  of  small  dissipation, 
whereas  this  is  not  the  case  for  a  generic  method.  This  is  an  important  consideration 
when  studying  near-ideal  situations  where  decoherence  is  small,  but  may  be  the  primary 
point  of  interest  or  limiting  factor  (as  in  studies  of  quantum  computation  in  real 
systems,  for  example).  Examples  of  such  systems  include  (i)  those  where  the  time 
scale  of  interest  is  very  long  compared  to  the  natural  scales  of  the  problem,  (ii)  systems 
with  sensitive  dependence  upon  initial  conditions  and  (iii)  systems  where  dissipation 
and  decoherence  are  small,  but  of  significant  importance.  Further  areas  of  future 
investigation  include  application  of  these  methods  to  such  systems  as  (i)  atomic  clocks 
of  various  types  [21],  and  (ii)  numerical  studies  of  decoherence  in  quantum  computers. 
This  latter  application  may  include  studies  of  decoherence-free  subspaces,  for  which  the 
structure  of  the  dissipation  operator  is  of  significant  importance. 
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LABORATORY  OPERATIONS 


The  Aerospace  Corporation  functions  as  an  “architect-engineer”  for  national  security  programs, 
specializing  in  advanced  military  space  systems.  The  Corporation’s  Laboratory  Operations  supports 
the  effective  and  timely  development  and  operation  of  national  security  systems  through  scientific 
research  and  the  application  of  advanced  technology.  Vital  to  the  success  of  the  Corporation  is  the 
technical  staff’s  wide-ranging  expertise  and  its  ability  to  stay  abreast  of  new  technological 
developments  and  program  support  issues  associated  with  rapidly  evolving  space  systems. 
Contributing  capabilities  are  provided  by  these  individual  organizations: 

Electronics  and  Photonics  Laboratory:  Microelectronics,  VLSI  reliability,  failure  analy¬ 
sis,  solid-state  device  physics,  compound  semiconductors,  radiation  effects,  infrared  and 
CCD  detector  devices,  data  storage  and  display  technologies;  lasers  and  electro-optics,  solid 
state  laser  design,  micro-optics,  optical  communications,  and  fiber  optic  sensors;  atomic 
frequency  standards,  applied  laser  spectroscopy,  laser  chemistry,  atmospheric  propagation 
and  beam  control,  LIDAR/LADAR  remote  sensing;  solar  cell  and  array  testing  and  evalua¬ 
tion,  battery  electrochemistry,  battery  testing  and  evaluation. 

Space  Materials  Laboratory:  Evaluation  and  characterizations  of  new  materials  and 
processing  techniques:  metals,  alloys,  ceramics,  polymers,  thin  films,  and  composites; 
development  of  advanced  deposition  processes;  nondestructive  evaluation,  component  fail¬ 
ure  analysis  and  reliability;  structural  mechanics,  fracture  mechanics,  and  stress  corrosion; 
analysis  and  evaluation  of  materials  at  cryogenic  and  elevated  temperatures;  launch  vehicle 
fluid  mechanics,  heat  transfer  and  flight  dynamics;  aerothermodynamics;  chemical  and 
electric  propulsion;  environmental  chemistry;  combustion  processes;  space  environment 
effects  on  materials,  hardening  and  vulnerability  assessment;  contamination,  thermal  and 
structural  control;  lubrication  and  surface  phenomena. 

Space  Science  Applications  Laboratory:  Magnetospheric,  auroral  and  cosmic  ray 
physics,  wave-particle  interactions,  magnetospheric  plasma  waves;  atmospheric  and 
ionospheric  physics,  density  and  composition  of  the  upper  atmosphere,  remote  sensing 
using  atmospheric  radiation;  solar  physics,  infrared  astronomy,  infrared  signature  analysis; 
infrared  surveillance,  imaging,  remote  sensing,  and  hyperspectral  imaging;  effects  of  solar 
activity,  magnetic  storms  and  nuclear  explosions  on  the  Earth's  atmosphere,  ionosphere  and 
magnetosphere;  effects  of  electromagnetic  and  particulate  radiations  on  space  systems; 
space  instrumentation,  design  fabrication  and  test;  environmental  chemistry,  trace  detection; 
atmospheric  chemical  reactions,  atmospheric  optics,  light  scattering,  state-specific  chemical 
reactions  and  radiative  signatures  of  missile  plumes. 

Center  for  Microtechnology:  Microelectromechanical  systems  (MEMS)  for  space 
applications;  assessment  of  microtechnology  space  applications;  laser  micromachining; 
laser-surface  physical  and  chemical  interactions;  micropropulsion;  micro-  and  nanosatel¬ 
lite  mission  analysis;  intelligent  microinstruments  for  monitoring  space  and  launch  sys¬ 
tem  environments. 

Office  of  Spectral  Applications:  Multispectral  and  hyperspectral  sensor  development; 
data  analysis  and  algorithm  development;  applications  of  multispectral  and  hyperspectral 
imagery  to  defense,  civil  space,  commercial,  and  environmental  missions. 
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